The differential energy distribution of the universal density profile 

of dark halo 
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ABSTRACT 

We study the differential energy distribution of dark matter halos, carrying out cosmological 
TV-body simulation. We give an analytical formula of the differential energy distribution of 
dark matter in the halos obtained by the numerical simulation. From the analytical formula we 
reconstruct the density profile described by the Navarro, Frenk, & White (NEW) profile. The 
NFW profile is consistent with the analytical formula of our fractional mass distribution. We 
find that a parameter in our analytical formula of differential energy distribution which is related 
with the slope of inner cusp of dark halo. We obtain the distribution function for the NFW 
profile which has sharp cut off at the high binding energy. We discuss physical reason of form of 
the analytical formula. 

Subject headings: Cosmology: dark matter; Galaxies: Formation, Halos, clusters 



1. Introduction 

It is interesting to understand how the den- 
sity profiles of galaxies and clusters of galaxies 
have formed. Navarro, Frenk and White (1995, 
1996, 1997; NFW) have shown in their N-body 
simulations of Cold Dark Matter (CDM) in the 
standard biased CDM and four power law spec- 
tra with indices n = 0, —0.5, and —1, open CDM 
(f2o = 0.1) with power-law spectra {n — Q and 
— 1), and ACDM cosmology that density profiles 
of dark halos have an universal profile described 
as 
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Several investigators have shown that the for- 
mula provides a good fit to their numerical results 
(Cole and Lacey 1996, Tormen, Bouchet & White 
1997, Huss, Jain & Steinmmetz 1999, Thomas 
et al. 1998). However, some other simulations 
(Fukushige & Makino 1997; Moore et al. 1998; 
Okamoto & Habe 1999, 2000) indicate that their 
density profiles have steeper inner cusp than the 



NFW profile. Jing (2000) gives his numerical re- 
sults that steeper cusps of density profile are found 
in recent mergers and in dark halos with substruc- 
tures, by a large set of high resolution cosmological 
simulations. 

Subramanian, Cen and Ostriker (2000) dis- 
cussed the general theoretical grounds of the den- 
sity profiles. They discussed that there is a pos- 
sible connection between slope of inner region of 
dark halo and the formation epoch and proposed 
that there is possible relation to cosmological pa- 
rameters. 

However, physical reason of the NFW profile 
have not been explored, although it is very im- 
portant to understand formation of galaxies and 
clusters of galaxies. 

We study the differential energy distribution of 
the NFW profile, dM/de, where e is the binding 
energy. The differential energy distribution was 
studied to understand virialization process of TV- 
body system (van Albada 1982). Binney (1982) 
shows that dM/de of the elliptical galaxies is ap- 
proximated by the Boltzmann factor, if surface 
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brightness of elliptical galaxies obeys the de Vau- 
coureurs' r^l^ law and their mass-to-light ratio is 
constant. It is interesting to study phase space 
structure of virialized objects with the NFW pro- 
file, since the differential energy distribution may 
give insight of the relaxation process (Binney & 
Tremain 1987). 

In this paper we calculate formation of clusters 
of galaxies by ./V-body simulation and we study 
the differential energy distribution of the clusters 
in our numerical results. We find an analytical 
formula of the fractional mass distribution that 
is defined by the differential distribution divided 
by mass of the object, is fitted by an analytical 
formula. We show that the NFW profile is consis- 
tent with our analytical formula of the fractional 
mass distribution. Using the iteration method, we 
construct the density profile from the analytical 
formula to show how the slope of the cusp of the 
density changes with parameters in our analytical 
formula. 

In the next section, we illustrate our numerical 
method. In the section 3, we present our numerical 
results. We show the iteration method and give 
results by this method in the section 4, and we 
summarize and discuss our results in the section 
5. 

2. Numerical simulation 

We use the Hofman and Ribak's (Hofman and 
Ribak 1991) procedure to set initial conditions in 
order to have a massive dark halo near the cen- 
tre of a simulation box. The cosmology is SCDM 
model (e.g. Davis et al. 1985) (Q = 1, as = 0.67 , 
Ho = lOO/i km s"^ Mpc'Sand h = 0.5). 

Numerical simulations are carried out using 
GRAPESPH code. GRAPE is a special purpose 
hardware to calculate gravitation between N-body 
particles (Sugimoto et al. 1990). We combined 
Smoothed Particle Hydrodynamics (SPH) (Mon- 
aghan, 1992) with GRAPE. We select massive ha- 
los of which mass is as large as that of cluster of 
galaxies and calculate their density profile and the 
fractional mass distribution. 

Mass of a CDM particle and a SPH particle are 
5.89 X 1O"M0 and 3.10 x IO^OM©, respectively. 
Gravitational softening length is 100 kpc. Both 
number of CDM and SPH particles are 29855, re- 
spectively. Size of the simulation box is 80 Mpc. 



3. Numerical results 
3.1. Density distribution 

Before we show our numerical results, we sum- 
marize the NFW profile. NFW proposed that the 
profile of the dark halo of cosmological object in 
their numerical results as 



where = r^oo/c, c is a dimensionless parame- 
ter and Per = (3i?o ) / (SttG) is the critical density 
of the universe . Since the mass of the halo is 
-^(''200) = 200 X 47r/3/9CT''200' there is a relation 
between 8c and c as 

200 

'^'="^[ln(l + c)-c/(l + c)]' 

where M(r 200) is the mass of which averaged den- 
sity inside r2oo is 200 times Pcr- 

This density profile has a cusp in the inner re- 
gion, p{r) oc r~^. p{r) oc in the outer region. 
NFW97 showed that the parameter c decreases 
with halo mass. 

Table 1 is physical values of our simulated clus- 
ters. The units of mass, and X-ray temperature 
arc IO^^Mq and 10*^ K, respectively. 

Figure 1 shows a density profile of our simu- 
lated typical rich cluster, CLc, and the NFW pro- 
file with c = 4.4 which fits well the numerical re- 
sult. Density profiles of dark halos obtained by us 
agree well with the NFW profile in the range from 
the gravitational softening length to r2oo- 

3.2. The differential energy distribution 

We introduce the differential energy distribu- 
tion, dM/de which gives the mass of dark matter 
in the dark halo with binding energy between e 
and e + de, where e is the specific binding energy, 

£ = *(r)-^^;^ (4) 

and a relative potential, = — $ 4- $o- ^ is grav- 
itational potential and we choose $0 to be such 
that a distribution function, /, is / > for e > 
and / = for e < 0. In our analysis, s is normal- 
ized by GM2oo/''20o- And we also introduce the 



2 



fractional mass distribution as the differential en- 
ergy distribution divided by the total mass of the 
dark halo, N{e) = dM/de/M. 

Figure 2 shows the fractional mass distribution 
of CLc. In figure 2, we also show N{e) given by 

9/(1-9) 



where 



N{e) = No 



(5) 



with q = 0.667 and Eq = 1.47. Figure 2 shows that 
equation (5) agrees well with our numerical results 
in the range of 0.5 < e < 4. There is cut off near 
£ ~ 4. We find the fractional distribution N{£) 
can be fitted by following formula, for q ~ 0.6 — 0.7 
and Eq — 1-2 — I -50 M200 / 1^200 for rich clusters in 
our numerical results. 

In table 2 we give A^O) and so of our numerical 
results for rich clusters. 

We have shown from our N- body simulation 
that the fractional mass distribution is also well 
fitted by the equation (5). However, near s = 
4 in figure 2 there are small number of iV-body 
particles. We should confirm consistency between 
A''(e) given by equation (5) and the NFW profile. 
We give a fractional mass distribution from the 
NFW profile as follows for the comparison. 

We assume that phase-space distribution func- 
tion /(x, v) depends e. At a radius r, velocity of 
a dark matter particle of the binding energy, e, is 
V = a/2(^' — e). The density profile may be given 
as follows (Binney and Tremaine 1987) 



p(r)=47r / ' f{e)[2{^-et'^de, 



(6) 



where Vg is the edge of the dark halo. From this 
equation, we may give /(e) as 



1 d r 



dp/d^ 



x/Stt^ ds y,^,^ [e - *]V2 

where Smin = ^'(rg). 

Equation (6) gives mass M as 
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From equation (8), the differential energy distri- 
bution is 

dM{e) 



g{e) = 167r2 f "^"^[2(4- - e)Y'^r''dr, (10) 
Ja 

and rmi^) is maximum radius that can reached by 

a particle of the binding energy e. 

If we assume the density profile is the NFW 
profile, we get dM/de from the equations (7), (9), 
and (10) for the NFW profile. 

Figure 3 shows the fractional mass distribution 
of NFW obtained in this way and N{e) given by 
equation (5). We show that N{e) given by equa- 
tion (5) is consistent with the NFW profile. We 
note that the cut off at high binding energy seen 
in figure 2 is not artifact due to limitation of our 
numerical resolutions. 

Figure 4 shows the distribution function of the 
NFW profile obtained by the above method. Since 
lower binding energy particles evaporate from the 
dark halo, their fraction become small. This form 
of function is similar to the distribution function 
of King model, //f (e) oc {e^'"^ - 1) (Binney 1982). 
On the other hand, there is a peak and sharp cut 
off at the high binding energy. This part corre- 
sponds to the central cusp. Although the distri- 
bution function of NFW profile have cut off at 
the high energy, the distribution function of King 
model does not have such cut off. This is an im- 
portant difference between them. 

4. The fractional mass distribution, den- 
sity profile, and the distribution func- 
tion 

We find that the NFW profile satisfies the frac- 
tional mass distribution given by equation (5) with 
g ~ 0.6 — 0.7. We study how the density profile 
changes when we change the parameter q and £0 
in equation (5). In this study, we use an iteration 
method as shown in the next subsections. 

4.1. The iteration method 

Binney (1982, and see also Binney and Tremaine 
1987) studied the phase space structure of galaxies 
of which surface brightness is the de Vaucouleurs' 
r^/'^ law. We apply his method to our study of 
the phase space structure of dark halo with the 
NFW profile. We obtain the density profile and 
the phase space distribution which are consistent 



de 



(9) 
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with equation (5), using the iteration method. 
The procedure is as follows. 

We assume the dark halo is spherically sym- 
metric. From equation (5), we assume dM{e)/de 
given by equation, 



dM{e) 
de 



1 



(!-,)(- 



9/(1-9) 



(11) 



where Mq = 1.0. Wc use a trial fimction of a den- 
sity profile, to calculate a gravitational potential 
in 50 (e), 



1 + 



3/2 ' 



(12) 



where Xc = TcItiqq = 0.1. This trial function is 
similar to the /3-model of X-ray surface brightness 
of clusters of galaxies. We obtain relative poten- 
tial from this trial density function. The relative 
potential is written as 



*(a;) = 47rG 



+ 



\ r 2 

— / po{x')x' dx' 
X Jo 

/ po{x')x'dx' 

J X 



+ (13) 



where x = r/r-ioo, and Xo = 10 is assumed, 
is chosen to be / > for e > and / = for 
e < 0. We obtain gois) by equation (10) for pQ. 
Prom equation (9), 



foie) = 



dM{e)/de 
50 (e) 



(14) 



Next, we calculate new density as, 



Pi+i{x) = (1 - a)pi{x) 

+ax4Tr [ fi{eW2{^i-e)de, (15) 

where i is an iterative index and a is an arbitrary 
constant of convergence. Here we assume a = 
0.5. Next we again follow the step from equation 
(13) to equation (15) by using pi instead of po 
and we obtain p2. We continue until pi becomes 
converge. In this way, we obtain p and /(e) which 
are consistent with equation (11). 



4.2. The density profile and the distribu- 
tion function 

Using the Binney's iteration method, we recon- 
struct mass density profile from equation (11). 

Wc show our result for q = 0.67 and e = 1.4 
in figure 5. A solid line is given by the itera- 
tion method and a dashed line is the NFW profile. 
We confirm that dM{e)/de characterizes well the 
NFW profile. 

In figure 6, we show density profiles for different 
q but e = 1.4. Smaller q (e.g. q = 0.5) results in 

shallower core in the inner region. On the other 
hand, larger q {q > 0.67) makes a cusp steeper 
than the NFW profile, density profile approaches 

p oc r^^ in the inner part, for q — > 1. 

For various values of Eq, the density profiles are 
similar to the NFW for q = 0.6 as shown in figure 
7. Absolute value of the density depends on Eq. 
Therefore, the slope of the cusp depends on only 
q, not Eo- 

Figure 8 shows the distribution function, /, ob- 
tained by the iteration method for various values 
of q . These curves show the same dependence on 
e in < e < 1 . Peak values of / are different each 
other. We also show the Boltzmannian distribu- 
tion for comparison in figure 8. For large q, peak 
value of /(e) and the maximum binding energy 
of the distribution become large. We have shown 
that the density profile with large q have the steep 
cusp. Therefore, sharp peak of /(e) corresponds 
to the steep cusp. 

Figure 9 shows / obtained by the iteration 
method for same q but various eo. These curves do 
not show the same dependence on binding energy. 
Height of peaks of these curves is constant. 

5. Summary and Discussion 

Wc analyze the universal density profile of dark 
halo proposed by NFW and its differential energy 
distribution. Our main results are summarized as 
follows. 

1. We study the fractional mass function A''(e) 
for dark halo obtained by our numerical sim- 
ulation and find its analytical formula which 
is the equation (5). 

2. We show that the NFW profile is given by 
the equation (11). 
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3. We show that the slope of the cusp in the 
density profile changes with a value of the 
parameter q in the analytical formula. 

We can regard that N{e) shows the statistical 
property of the NFW profile. If the NFW profile is 
universal, q = 0.6 — 0.7 in equation (5). Different 
q makes slope of a cusp different. Since q plays 
an important role, we should make clear what 
physical process determines q. Recent high resolu- 
tion numerical simulation (Okamoto & Habe 1999, 
2000) shows the steeper cusp, p oc r~^'^, than 
the the NFW profile. This profile corresponds to 
q = 0.75 — 0.8. Isothermal profile, p ex r^^, corre- 
sponds to g = 1. 

We study /(e) for the NFW profile. The for- 
mula of this is not isothermal one nor the King 
formula /k oc e^/-^' - 1. /(e) for the NFW profile 
have the energy cut off at the high end of e. We 
should study the reason why /(e) has such a form. 
Lynden-Bell (1967) studied distribution function 
/(e) of a virializcd system. Maximizing the Boltz- 
mann entropy of the system, resulting distribution 
is isothermal profile, p oc r~^. In this case the sys- 
tem has infinite extend and infinite mass. This is 
not realistic for astronomical objects. Cosmologi- 
cal simulations have shown that galaxies and clus- 
ters of galaxies formed in these simulations have 
more rapid radial decline than isothermal in the 
outer part. 

We note that the form of the distribution func- 
tion relates with the cusp profile. For q < 1, 
there are the maximum binding energy Cmax and 
the sharp peak in /(e) at emax- Since there is 
£max, phase space volume occupied by dark mat- 
ter is limited than the isothermal profile. This 
may be the reason of the cusp profile different from 
p oc r~^. For g < 1, the distribution function /(e) 
changes with q values as shown in figure 8. 

We note that the form of equations (5) and (11) 
are similar to the Tsallis' escort distribution, 



fractral space time) the following generalized en- 
tropy has been proposed : 



P{E,T') 
P{0,T') 



9/(1-9) 



(16) 



where T' is temperature parameter and q is en- 
tropic index (TsaUis, Mendes, & Plastino, 1998). 
Tsallis' non-extensive generalized statistics (Tsal- 
lis 1988) is paid attention in the area of statistics 
of a multi-fractal system. In non-extensive system 
(long-range microscopy memory, long range forces. 




where fc is a positive constant. Optimization of Sg 
yields, for the canonical ensemble, 

Pi = Z;'[l-{l-q)e,/T']'/^'-"\ (18) 
Z,^^[l-(l-g)e,/r']i/(i-«) (19) 

and, when q — > 1, the Boltzmann-Gibbs result is 
recovered. In this statistics an expected value of 
any physical variable is given by the Tsallis' escort 
distribution : 



E 



3 ^3 



(20) 



where {Ai} are the eigenvalues of an arbitrary 
observable A. The escort distribution, Pi = 
pf/^jP'j, is similar to the form of equation (5) 
and (11). 

It is expected that long range interaction makes 
the system non-extensive. Our case may be the 
one of this cases. Our results suggests that dif- 
ferential energy distribution of collisionless parti- 
cles is described by the escort distribution of the 
Tsallis statistics. This may indicate that energetic 
process of collisionless particles must be stochastic 
process. 

In the Tsallis' escort distribution function, 
there is the maximum value of e for pi > Q for 
< g < 1. This case is called superextensive. 
The shallower cusp profile in the the NFW pro- 
file shows the superextensive property of the dark 
matter distribution. 

We show that the NFW profile corresponds to 
g = 0.6 — 0.7. Lavagno et al. (1998) have recently 
shown that fraction of peculiar velocity of cluster 
of galaxies (Bahcall & Oh 1996) is well explained 
by the Tsallis escort integral, 

X - (1 - 'z)(«/«o)')«/^'-'^^dw 

~ jy-{l-(l-q){v/vor)'^/('-'^->dv- 

They obtained q = 0.23 to fit the fraction of pe- 
culiar velocity of cluster of galaxies and is smaller 
than in our case. There is a conjecture that q 
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of system approaches unity when the system pro- 
ceeds relaxation (Tsallis 1999). Our results in 
which = 0.6 — 0.7 are consistent with this con- 
jecture, since dark matter distribution in a clus- 
ter of galaxies is more relaxed system than the 
large scale motion of clusters of galaxies. The 
value of q must bo related with the degree of relax- 
ation of coUisionless particles, g = 1 corresponds 
to isothermal distribution in the Tsallis statistics. 
The gravitational system like clusters of galaxies 
with the universal profile may be non-extensive 
because they formed recently. 

We should study the reason why dark halo has 
the value of g = 0.6 — 0.7 in the hierarchical clus- 
tering scenarios. It is interesting to study the dif- 
ferential energy distribution of self-gravitational 
system formed in a circumstance without hierar- 
chical clustering to make clear mechanism what 
determines q of the gravitational system. 

We would like to thank Sumiyoshi Abe, Masayuki 
Fujimoto, and Seiichi Yachi for helpful comments 
and discussion. 
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Fig. 1. — Density profile of one of simulated clus- 
ters, CLc, and the NFW profile with c = 4.4. Solid 
line is our numerical result and dashed line is the 
NFW profile. Radial distance is normalized by 
r2oo- Our simulated clusters agree well with the 
NFW formula between the gravitational softening 
length to r2oo- 
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Fig. 3. — Fractional mass distribution of dark 
halo. Solid line is caliculated by equations (7), 
(9), and (10) and dashed line is given by equation 
(5) with q = 0.66 and sq = 1.4. 
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Fig. 2. Fractional mass distribution of our simu- 
lated rich cluster, CLc. Solid line is our numerical 
result and dashed line is given by equation (5) with 
q = 0.667 and eo = 1-47. 
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Fig. 4.— f{e) of the NFW profile. Solid line is 
/(e) and dashed line is a fitted exponential func- 
tion. 
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Fig. 5. Density profile constructed by the iter- 
ation method, assuming the equation (11). Sohd 
Hne is our result for q = 0.67 and so = 1.4 and 
dashed line is the NFW profile of c = 6.59. Radial 
distance is normalized by r2oo- 
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Fig. 7.- Density profiles for various e. Dotted, 
short-dashed, dashed, and solid lines are for Eq = 
0.5, 1.0, 1.5, and 2.0, respectively. Density profiles 
are self-similar. However, their normalizations are 
different. 
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Fig. 6. — Density profiles for various q. Dot- 
dashed, dotted, short-dashed, dashed, and solid 
lines are for q = 0.25,0.5,0.67,0.75, and 1.0, re- 
spectively. Density profile with smaller q has a flat 
core. On the other hand, one approaches p oc 
in r < 0.1 for q — > 1. 
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Fig. 8. — We show /(e) for various q ; Dot-dashed, 
dotted, short-dashed, dashed, and solid lines are 
for q = 0.25,0.5,0.67,0.75, and 1.0, respectively. 
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Fig. 9. — We show f{s) for various sq; Dotted, 

short-dashed, dashed, and sohd hnes are for eo = 
0.5, 1.0, 1.5, and 2.0, respectively. 



Table 1: The physical vahies of our sinmlatcxl clus- 
tors 







Mass 
(lOi^Mo) 


T 

(10^ K) 


c 


CLa 




3.05 


0.953 


2.50 


CLb 




2.80 


0.890 


2.72 


CLc {z 


= 0) 


3.36 


1.29 


4.41 


CLc {z 


= 0.25) 


1.60 


0.796 


2.31 


CLc [z 


= 0.5) 


0.831 


0.540 


2.67 


CLd 




3.71 


1.57 


6.34 


CLe 




2.82 


1.19 


4.65 


CLf 




2.31 


0.654 


2.75 



Note. — Each value is calculated at redshift z = 0, except for CLc. 



Table 2: iVo, q, and eq of rich clusters in our nu- 
merical results 





No 


9 


£o 


CLa 


0.309 


0.612 


1.27 


CLb 


0.253 


0.680 


1.21 


CLc (z = 0) 


0.190 


0.667 


1.47 


CLc (z ^ 0.25) 


0.230 


0.606 


1.05 


CLc (z = 0.5) 


0.305 


0.600 


0.672 


CLd 


0.276 


0.699 


1.26 


CLe 


0.117 


0.503 


2.03 


CLf 


0.258 


0.580 


1.41 



Note. — Each value is calculated at redshift 2 = 0, except for CLc. 
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